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Abstract 

Deformations of liquid interfaces by the optical radiation pressure of a focused laser 
wave were generally expected to display similar behavior, whatever the direction of 
propagation of the incident beam. 

Recent experiments showed that the invariance of interface deformations with re- 
spect to the direction of propagation of the incident wave is broken at high laser 
intensities. In the case of a beam propagating from the liquid of smaller refractive 
index to that of larger one, the interface remains stable, forming a nipple-like shape, 
while for the opposite direction of propagation, an instability occurs, leading to a 
long needle-like deformation emitting micro-droplets. While an analytical model suc- 
cessfully predicts the equilibrium shape of weakly deformed interface, very few work 
has been accomplished in the regime of large interface deformations. In this work, 
we use the Boundary Integral Element Method (BIEM) to compute the evolution of 
the shape of a fluid-fluid interface under the effect of a continuous laser wave, and 
we compare our numerical simulations to experimental data in the regime of large 
deformations for both upward and downward beam propagation. We confirm the 
invariance breakdown observed experimentally and find good agreement between 
predicted and experimental interface hump heights below the instability threshold. 
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1 Introduction 



The deformation of liquid-liquid interfaces by the optical radiation pressure 
has received increasing attention in the past few years as many practical ap- 
plications of laser-induced surface deformation are now under development. 
Among others, we can cite interfacial characteristics measurements such as 
viscosity jl| or surface tension as well as fluid membranes manipulation 
with optical tweezers 0; 0|. Effects of the radiation pressure have been also 
recognized as an appealing non-intrusive tool for local manipulation of liquid 
or soft materials giving birth to many applications in biotechnologies 

Historically, the deformation of fluid-fluid interfaces resulting from the radia- 
tion pressure induced by an impinging focused laser beam was flrst identifled 
in the early 70's by Ashkin and Dziedzic Later, Zhang & Chang il- 
luminated a 50 /im radius water drop with a 100 — 200 mJ laser pulse and 
showed strong distortion of the droplet surface at its front and rear regions. At 
low pulse energy, the droplet interface exhibited oscillations that damped out 
on time scales of tens of microseconds, whereas at high energy a disruption 
generating a jet of micro-droplets was observed at the rear part of the drop. 
The drop distortion was theoretically studied by Lai et al. Q and later by 
Brevik et al. Based on a linear wave theory adapted to low energy pulses, 
their analysis predicted drop oscillations very similar to those observed in the 
experiments. However, drop deformations under higher energy pulses were not 
modelled as their linear model can no longer be used to describe the droplet 
shapes in the regime close to disruption. 

Thus, current existing theoretical descriptions of optically induced flow and 
surface deformations are restricted to small amplitude deformations. However, 
recent experiments on very soft interfaces have evidenced several regimes, 
ranging from the so-called classical linear regime for small beam intensity, 
in which the height of the deformation linearly depends on the radiation to 
Laplace pressure ratio referred to in the following as ^ to nonlinear regimes 
with a possible interface breakup for even larger beam intensities 14 



In these recent experiments, a continuous Ar"*" laser wave (wavelength in vac- 
uum Ao = 514.5 nm) of waist cuq ~ 3 to 15 /im was used to bend the interface 
between two liquid phases in coexistence close to their liquid-liquid critical 
point. Two reasons motivated this choice. As the separated phases of these 
near-critical mixtures have very low surface tension (7 ~ 10~^ N/m), a laser 
beam of moderate power P ~ 1 W becomes sufficient to induce large inter- 
face deformations of typical size ~ 10 — 100 /im. Moreover, near-criticality 
raises universality concepts demonstrating the generality of the purpose. As 
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predicted from the photon momentum balance at the interface, in the hnear 
regime the radiation pressure induces the same deformations for upward or 
downward laser propagation 



11; 13 



However, in the nonlinear regime the 
invariance of the interface deformation with respect to the direction of propa- 
gation of the wave breaks down. Stable nipple-like deformations were observed 
in the case of a propagation from the less refractive fluid (marked as (1)) to 
the more refractive one (marked as (2)) {upward propagation) while in the 
other case where the laser wave propagates from the more refractive fluid to 
the less refractive one (downward propagation) , the interface becomes unstable, 
leading to the formation of a needle-like deformation emitting micro-droplets 
(3 (see Fig. 1). 



(a) 





Fig. 1. Interface deformations induced at (T — Tq) = 3.5K by a laser beam of waist 
(jQ = 5.3fim. (a) Laser propagating upward from the less refractive fluid to the more 
refractive fluid as indicated by the white arrow. P increases from top to bottom and 
is successively equal to 120, 240, 360, 390 and 720 mW. (b) Downward direction of 
propagation. P= 124, 248 and 372 mW from top to bottom; the bottommost picture 
(405 mW) shows the destabilization of the interface leading to the formation of a 
stationary jet similar to that illustrated in (c). 

The aim of the present work is to investigate whether the differences in in- 
terface deformation with respect to the direction of propagation of the wave 
(called hereafter invariance breakdown) is numerically predictable and if the 
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predicted interface shapes and heights agree with experiments. Our paper 
is structured as follows. In Section 2, we briefly describe the experimental 
setup, sample properties, and the physical model for laser-induced interface 
deformations. This model assumes an axisymmetric Stokes flow in each liquid 
and a boundary condition at the interface describing the competition between 
viscosity, optical radiation pressure, capillarity and gravity effects. A brief de- 
scription of the Boundary Integral Element Method (BIEM) is presented in 
Section 3, emphasizing the advantages of its application to interfacial flows. In 
Section 4, comparisons between numerical and experimental results are shown 
and discussed for both the linear/nonlinear regimes of deformation and both 
directions of propagation, illustrating by the way the efficiency of our approach 
to investigate the subtle coupling between the effects of light and flow. 



2 Experimental configuration and physical model 



Exhaustive experimental details on the conflguration and protocol used here 



11; 13 



were reported earlier 



2. 1 Experiments 



In figure 2 we have represented a picture of a typical liquid-liquid interface 
deformation induced by a laser beam propagating upward, together with the 
notations used throughout this work. Cylindrical coordinates (er, Bz, Bq) with 
their origin O located at the intersection of the beam axis with the initial fiat 
interface are chosen for this study and are shown in Fig. 2. A point x is thus 
marked by the space coordinates (r, 

Incidence and transmission angles of light are respectively denoted by Qi and 

The two-phase liquid sample is enclosed in a fused-quartz cell of optical path 
length / = 2 mm. The bending of the liquid-liquid meniscus is driven by a 
linearly polarized continuous Ar^ laser (wavelength in vacuum Aq = 514.5 nm) 
in the TEMqo gaussian mode. The beam is weakly focused on the interface by 
a lOx microscope objective to ensure a cylindrical distribution of the intensity 
near the meniscus. Thus, in the vicinity of the liquid-liquid interface, the laser 
beam intensity I has the following distribution: 



/(r, z)^/(r) = ^exp (-2 f— ) 
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Fig. 2. Liquid-liquid interface deformation by a focused laser wave propagating up- 
ward. T — Tc = 3.5K; ujq = 5.3^m ; P = 240mW. The laser light scattered towards 
the camera has been filtered. See text for notations. 

where P is the laser beam power. The beam waist ujq can be adjusted from 3 
to 15 yUm. Depending on the setup configuration, the laser beam propagates 
vertically either upward from fluid (1) to fluid (2) or downward from fluid (2) 
to fluid (1). 



2.2 Two-phase sample properties 



The investigated fluid-fluid interface is obtained according to the following 
procedure. Using a quaternary liquid mixture made of toluene, sodium do- 
decyl sulfate, n-butanol and water, we prepare a water-in-oil micellar phase 
of microemulsion whose composition is such that at room temperature it is 
close to a critical consolute line. Close to the liquid-liquid critical temperature, 
Tc = 308 K, the critical thermodynamic behavior of the mixture belongs to 
the universality class (d=3 n=l) of the Ising model For a temperature 
T > Tc it separates in two micellar phases of different concentrations ipi and 
ip2- The use of a near-critical two-phase fluid is motivated by the fact that 
significant interface deformations by electromagnetic radiation, without non- 



linear propagation effects [l5| or disturbing thermal coupling, require weak 
interfacial tension and buoyancy. With our fiuids, both effects are fulfilled 
since they vanish when the critical point is neared. Moreover, the interfacial 
tension 7 of phase-separated near-critical supra molecular fiuids is intrinsi- 
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cally smaller than that of near-critical pure fluids, thus enhancing even more 
interface deformations. 



The mixture is thermally controlled at a temperature T above Tc- Since the 
density (index of refraction) of water is larger (smaller) than that of toluene, 
the micellar phase of largest concentration ipi, density pi and smallest refrac- 
tive index Ni is located below that of lower micellar concentration (y92, density 
P2 and of refractive index A''2. Given the temperature range investigated in 
the reported experiments T — Tc = 2 to 15 K, thermophysical properties of 
the two-phase sample can be satisfactorily evaluated using usual asymptotic 
scaling laws of near-critical phenomena. For interfacial tension: 



T-T, 



2v 



7 = 7o(^^) , (2) 



with 7o = 10~^ N/m and z/ = 0.63. Assuming that (i) a scaling law accurately 
describes the variations of Ay? = — ^2 in the investigated broad tempera- 
ture domain T — Tc and that (ii) the coexistence curve is symmetric versus 
the critical concentration yjc, the micellar concentration in each phase can be 
estimated by: 



(3) 
(4) 





Aip 




Alp 


with (pc = 


= 0.11 and: 


Aip -- 


= A^, (— 



(5) 



with (3 = 0.325. 



The value of the critical amplitude A{pc, can be estimated theoretically for 
this system with the main assumption that our mixture is binary Aipc 

sions 



18 



1.458 16; 17 



R'^ = 0.37 being a universal ratio in microemul- 
As micellar phases are in fact quaternary components fluids, their 
phase diagram presents some asymmetry in ip leading either to over or under- 
estimate the critical amplitude Ay^c depending on which side of the coexistence 
curve is chosen. The theoretical value of Aipc leads to an overestimation of the 
optical index contrast AN. Consequently a modified value for Aip^, 0.42, has 
been adopted in this study, based on a quantitative comparison between nu- 
merical and experimental interface steady hump heights in the linear regime 
of deformations. This modified value gives acceptable predictions for the con- 
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centrations (fi and (p2- In fact, Aipc should be understood as a free parameter 
in the estimation of phases properties. 



The density of each phase pi, i = 1,2 can be written as a function of ifi: 

Pi = Pmic^i + PcontO- - Vi)-, (6) 

where pmic = 1045 kg.m"'^ and pcont = 850 kg.m~^ are the densities of the 
micelles and continuous phases respectively. At T — Tc = 3.5K the estimated 
density contrast is Ap = pi — p2 =— 20kg /m^. 

As the average distance between two micelles is small compared to the wave- 
length of the laser wave, the mixture can be regarded as homogeneous from 
the electromagnetic point of view. Thus, the mean-field model for the relative 
dielectric permittivity ej of the mixture predicts 



21 



ei{(pi) = i^iemic + (1 - <^i)econi - —, -71 ^ V- (7) 



This relation is used to estimate iVj taking into account the fact that: 

e, = iV2, z = l,2. (8) 



along with e^c = 1-86 and econt = 2.14, e^ic and econt being the relative 
permittivity of the micelles and continuous phases respectively. 

At T - Tc = 3.5ir, the optical indices are A^i = 1.447 and N2 = 1.457. 

In addition, since concentrations are weak, we use Einstein's relation to esti- 
mate the dynamic viscosity pi of each phase: 



p2 = po[l- 2.5^) , (10) 



/ii = /io (1 + 2.5^) (9) 
with po = 1.269 Pa.s. 

Nevertheless, these values could be shifted by possible laser heating of the 
fluids. In order to estimate the resulting change in the physical properties such 
as interfacial tension 7 or viscosity pi due to temperature increase, we consider 
the steady diffusion equation in cylindrical coordinate with the absorbed laser 
intensity as a source term: 

V^T,(r) + ^/(r) = 0, (11) 



th 
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We assume stationary conditions since the thermal diffusion time scale is much 
smaller than the viscous one. T/(r) is the local increase of temperature due to 
local heating of the laser wave, ath — 3 10~^cm~^ is the thermal absorption 
and Ath = 1-28 10~^Wcm~^ is the thermal conductivity. Using a Fourier- 



Bessel transform to solve equation (iTTIl |19l : Il6l |. we find that the maximum 
increase in temperature is: 



Ti{r = 0) ~ ^^IniWOT) (12) 



where T = 1.781 is the Euler constant. Considering now equation ([2]), we can 
estimate the change of interfacial tension due to the increase in temperature: 

dj 1.26 , , 



At T - Tc = ?>.hK and for P = IW, we find that ^ ~ 1.3 IQ-'^Nm-^K-^ 

A7 

and T/(r = O) ~ O.li^ which leads to ~ 3.6%. Thermocapillary effects 

7 

can thus be confidently discarded. 

Considering now the viscous dependance on temperature, we use the following 



empirical law given for microemulsions [20 



/i(T) = [1.934 - 0.019(T - 273)]10"^ (14) 



Consequently, = -0.019 10"^ and thus for P = IW we find — ^ 2%. 
oT n 

This second estimation ensures negligible thermal effects. 



2. 3 Electromagnetic force and pressure 



As expressed in [2l||, the total electromagnetic force per unit volume exerted 
by the laser in each phase is given by: 



-^eoEi^VQ + ^eoV 



dpi 



+ 



e,; - 1 a 



dt 



:EixHi). 



(15) 



In this expression, E (respectively H) is the electric (respectively magnetic) 
field associated to the laser wave, eo is the permittivity of vacuum and c = 
3 lO^m/s is the celerity of light in vaccum. 
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The first term, whose jump across the hquid interface results in the usually 
called optical radiation pressure, is due to the change in photon momentum 
from one fluid to the other. This momentum change is due to the disconti- 
nuity of permittivities across the interface. As both fluids are assumed to be 
homogeneous, the optical radiation term cancels within each phase and only 
acts on the interface. 

The second term, referred to as the electrostrictive force, also undergoes a 
jump at the interface due to the difference in both the optical properties and 
the electric flelds E between the two liquids in contact. However, this force 
also acts as a bulk force within each phase because of the radial dependence of 
the electric held E. We will demonstrate, in adequation with previous theoret- 
ical investigations that the electrostriction does not contribute to the 
motion and shape of the interface as its bulk contribution is compensated by 
its surface one. The quantity which depends on the density of each phase 
through a nonlinear relationship can be deduced from equations ([6]) and (JT]). 

Finally, in equation ( fTSi l. the third term is called the Abraham term, and is 
undetectable at optical frequencies and thus cancels out for our purpose. 

In the following, we denote E"^ =< E"^ > the quadratic value of the electric 
held averaged over an optical period. 

Assuming incompressible fluids {pi is homogeneous throughout each phase i, 
i = 1,2), we can include contributions of the bulk forces (gravity and elec- 
trostriction) in the pressure field and define a pseudo-pressure Pi given by: 

Pi = Pio + PiQZ - ^eo f-^'Vi^ j ,^ = 1,2, (16) 



where piQ, z = 1,2, is the pressure in each phase and g = 9.81 m.s ^ is the 
acceleration of gravity. 



The electromagnetic stress tensor defined by Landau [2]J such that V.Tf™ = fg, 
can be expressed as: 

Xr = \eo (E^^P^^] I - l^oe^E^H + eo^EiEi, z = 1, 2. (17) 



In the case of a steady interface, the jump from fluid 1 to fluid 2 of the 
local pressure and of the electromagnetic stress are balanced by the Laplace 
pressure: 

(P20 - Pio)n + [Tr - Tr].n = ^K{r)n, (18) 
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where n is the unit vector directed from liquid 1 to liquid 2 and normal 

Id rz' 

to the interface and K,[r) = —- , is the curvature of the interface in 

cylindrical coordinates, z' = ^ is the local slope of the interface. 



Rewriting Eq. (iTSjl in terms of pseudo-pressures leads to: 



/ \ / \ ^ / T-i 2 ^^2 ^^ 2 \ 

{P2 - Pi) + (pi - P2)9Z + -eo{E2 P2t: pij^) 

2 dp2 opi 

+n.[Tr-Tr].n = 7'^(r). (19) 

Therefore, the bulk contribution of the electrostriction is balanced by its sur- 
face contribution and the interface shape does not depend on electrostriction. 
Consequently, the equilibrium equation of the interface can be written as: 

(P2 - Pi) + (Pi - P2)gz - ^eo{eiEi^ - es^g^) + eo(eiEiEi.n - eaEaEa.n).!! = jK{r).{20) 



At final equilibrium, only normal stress act on the interface, consequently 
there is no fiow within the phases and this equilibrium in each phase is char- 
acterized by: 



Vpi = 0,1 = 1,2. (21) 



In the experiments, the laser wave was linearly polarized so that the electric 
field E = Ee^ was perpendicular to the plane of observation of azimuthal 
coordinate a = 0. Therefore E is continuous across the interface in this plane. 

Conversely, in the plane defined by a = 7r/2, the electric field is within the 
plane of propagation (parallel polarization) and thus is no longer continuous 
across the interface. However, as shown in [l2j, the first term of the electro- 
magnetic force, which is responsible for the optical radiation pressure acting 
on the interface, is quasi-independent of beam polarization as long as refrac- 
tive indices are sufficiently close to each others (^ ~ 1). We can then assume 



as in a previous investigation [22| a circular polarization of the electric field E. 



Therefore, using the definition of the irradiance /(r) = ^NicEi^ in the case 
of a beam propagating upward, the radiation pressure can be written as: 

n"P(r, Oi, Ot) = cos0i(2iVi cosOi - T^^^iNi cosdi + N2 cosOt)), (22) 
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while in the case of a downward beam propagation it is given by: 

U'^'""''{r,e„et) = ^ cosei{2N2Cose, - T'^'""''{N2Cose, + Nicoset)).{23) 



Here, T"^ and T'^""'" are the Fresnel transmission coefficients in energy. In the 
case of a circular polarization, they are expressed as: 



rpup ^ j.down ^ ^2NiN2 COS 9, COS 9t] 

( 1 



(iVi COS Qi + N2 COS QtY {N2 cos Oi + Ni cos Ot 



i2 



(24) 



As a final remark, note that the expression of the electromagnetic pressure 
in the case of a downward propagation is valid until achieving the condition 
of total refiection 6i<6TR = arcsin(^^^. What occurs when light is totally 
reflected by a highly deformed interface is out of the scope of the present 
numerical study. 



2.4 Flow equations and boundary conditions 



We consider a cylindrical domain of radius R and height H as depicted in 
Figure 3. We choose to treat all the governing equations in a dimensionless 
form by using the laser waist ujq as the characteristic length scale and the 
viscous relaxation velocity m* = 7//i2 as the reference velocity associated to 
the characteristic timescale r* = /i2Co'o/7- The reference pressure is taken as 
p* = , i = 1,2. Since the weak absorption of light at laser frequency 
ensures negligible thermal effects, we consider all liquid properties (7, pi, pi, 
i = 1,2) constant in adequation with our estimation of the change in inter- 
facial tension and viscosity due to a local increase in temperature. Moreover, 
given the value of the Reynolds number of the flow P2U*ujq/ 112 ~ 10"^, iner- 
tia is negligible compared to viscous forces, allowing a quasi-steady creeping 
flow assumption. Thus, flow in each liquid is governed by Stokes and mass 
conservation equations respectively given by: 



- Vpi + V^Ui = (25) 
V.Ui = , 2 = 1,2. (26) 

Ui is the dimensionless velocity in fluid i and Pi is the dimensionless value 
of the pseudo-pressure deflned in Eq. (fT6ll . The hydrodynamic divergence free 
stress tensor Ti is: 

Ti = -pJ + (Vui+*Vui). (27) 
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According to Eq. (flGll . the normal stress jump across the interface Sj, which is 
compensated by gravity, optical radiation pressure and capillarity forces can 
be expressed in the following form: 

(ATi.n - T2.n).n = ^(r) - U{r) - Bo z. (28) 



In this equation, A = /ii//i2 is the viscosity ratio and Bo = (pi — p2)gujQ/'y 
is an optical Bond number which represents the ratio of gravity to capillary 
forces, n(r) is the dimensionless expression of the optical radiation pressure 
resulting from equations (l22l) and (l23ll . We denote n(r) = n(r)"^a;o/7 for the 
upward direction of propagation and n(r) = n(r)'^°"'"a;o/7 foi^ the downward 
direction. In order to quantify the effects of the laser wave on the interface, 
we define ^ as the electromagnetic to Laplace pressure ratio on beam axis at 
normal incidence: 

e = n(r = o,e, = o,e, = o) = I (29) 



A no-slip condition at the interface along with the fact that fluids are immis- 
cible implies continuity of the velocity u on the interface Sj: 

u = Ui = U2 for X e Sj. (30) 



Moreover, we assume the classical no-slip boundary condition on the lateral, 
upper and lower solid walls Sci and Sc2 (see figure 3). This leads to: 

Ui = for X e 5ci , « = 1, 2. (31) 



In order to reproduce the experimental configuration, R must be large enough 
compared to ujq to meet the hypothesis of an infinite extent in the horizontal 
direction. Several tests were performed with different values of the dimension- 
less radius /3 = ^ of the computational domain, and a value /3 = 70 was found 
to be large enough to satisfy this constraint. The interface motion is obtained 
using a Lagrangian approach. It consists in tracking each fluid particle on the 
interface in its Lagrangian motion according to: 

^ = u(x) for X e Si. (32) 



That is, the interface is advected along with the flow until equilibrium is 
reached for which velocities normal to the interface are zero, i.e. u(x).n = 
for t too- 
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3 Numerical method 



A brief description of the numerical algorithm is presented in this section. For 
more extensive details on the Boundary Integral Element Method (BIEM) ap- 
plied to two-phase axisymmetric flow, the reader may refer to the review by 
Tanzosh et al. on the solution of free surface flow problems using this tech- 



nique [23||. The BIEM reveals to be an excellent tool to solve interfacial flow 
problems with high resolution as reported in the analysis of flow involving 
electric and magnetic flelds [2j| or buoyancy [25|; l26| . 



25 
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Fig. 3. Typical configuration of the computational domain showing the two fluids 
and the liquid interface. The domain is assumed to be axisymmetric along e^. 



Because solutions to the Stokes problem can be formulated in terms of Green's 
functions [27||, we can rewrite the governing equations as a system of integral 
equations over the boundaries of the computational domain. Once boundary 
conditions on Sj, Sci and Sc2 are used, the two-phase Stokes problem can be 
written in the following compact form: 



1 + A 



u(x) = J \J.n{K{ry) - n(ry) - Boz{ry))dSy+ 



Si 



;i-A) J n.K.udSy + X J V.{Ti.n)dSy- J V.{T2.n)dSy. 

Si Sci Sc2 



(33) 



Here, U and K are Green kernels for velocity and stress respectively and are 



given by 



27 : 
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3 ,ddd, , , 

where d = x — y, y(rj,, z^) is the integration point. In Eq. (ISS]) . the first term 
in the right hand side describes the fiow contribution from interfacial tension, 
radiation pressure and gravity, whereas the second term accounts for shear 
rates contrast on the interface. This term vanishes when A = 1. The third and 
fourth terms account for shear occurring on (S'ci) and [Sci] as a result of the 
no-slip boundary condition. 

Velocities on the interface as well as stress over all the boundaries Sj, Sqx 
and Sc2 are determined by solving the discrete form of this equation using a 
numerical procedure. This procedure requires first the discretization of all the 
boundaries 5/, S'ci and Sc2- Due to integral formulation and axial symmetry, 
the problem is reduced to one dimension and only line boundaries, as repre- 
sented in Figure 3, need to be discretized. In this work, the mesh is made of 
constant boundary elements i.e. line segments with centered nodes. The fiuid- 
fiuid interface Sj is parameterized in terms of arc length and is approximated 
by local cubic splines, so that the curvature can be accurately computed. Dis- 
tribution and number of points are adapted to the shape of the interface, so 
that the concentration of elements is higher in regions where the variation of 
curvature of the interface is large. The number of mesh points is about 70 for 
a typical computation of a small interface deformation. The solid boundaries 
Sc\ and S'ci are meshed using about 40 uniformly distributed points. An in- 
crease in the mesh resolution for the interface and the solid boundaries do not 
show any significant change in the results. 



Azimuthal integration of Eq. (I33l) is performed analytically [28|; l29fl reducing 
Eq. (l33l) to a line integration which is finally performed using Gauss quadra- 
tures |30y. Elliptic integrals resulting from the azimuthal integration are eval- 
uated using power series expansions 

The motion of the interface is followed using the kinematic condition ( l32l ) 
which is discretized using an explicit first-order Euler time scheme. A typical 
computation begins with a fiat interface at rest. The laser beam is switched on 
at t = 0, and the interface deforms towards fiuid 1 of smallest refractive index. 

Computation stops when an equilibrium state is reached ( — ~^ ^'^^ 
time step is chosen to be about 20 times smaller than r*. 



4 Results and discussion 



In this section, we compare equilibrium hump heights and interface shapes ob- 
tained experimentally to their numerical predictions in both linear and non- 
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linear regimes and for both directions of propagation. We also present the 
dynamics of hump formation expected in the nonlinear regime for both up- 
ward and downward propagations. 

As mentioned before, few work were dedicated to the modelling of optical de- 



formations of liquid interfaces in the nonlinear regime. Hallanger et al. [22 
used a finite difference method to solve the equation describing the equilib- 
rium state of the interface for a laser wave propagating from the fluid of 



smallest refractive index (here fluid 1) and Wunenburger et al. [13|| used a 
simple model where gravity was neglected (Bo << 1) to predict the equilib- 
rium hump heights in the nonlinear regime. 




Fig. 4. Variation of the reduced hump height h = \z{r = 0)| versus the reduced 
pressure ratio ^ when light propagates upward and downward for = 5.3/im and 
T — Tc = 3.5K {Bo = 0.015). Comparison is made between experimental results 
(symbols) and the numerical resolution (lines). The bold dashed line is the prediction 
for the linear model. ^^^^ and ^^""^ respectively represent the experimental and 
numerical thresholds of instability in the case of downward propagation. 



In Figure 4, we have represented the variation of the dimensionless equilibrium 
interface hump height h as a function of the pressure ratio ^ for upward and 
downward propagations. Experiments (symbols) and numerical predictions 
(lines) are shown together. 
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4-1 Steady-state deformations, linear regime 



Previous experiments l32l | performed in the linear regime of deformations 
showed that, for ^ < 2.5, /i varies linearly with ^ for both directions of propa- 
gation (symbols) and that the linear variation of h{^) does not depend on the 
direction of propagation for N1/N2 — 1. A linear theory assuming Ni ~ N2 
was proposed [uJ; [l3| which quantitatively agree with these observations. In- 
deed, the linear theory gives: 

z{r) = -i y ^-^^Mkr)kdk, (36) 



where Jq is the zeroth-order Bessel J function. We deduce: 



c Bo „ . Bo . 

h{r = 0) = ^-e-E,{—), (37) 



— k 

where Ei{x) is the one- argument exponential function [Ei(x) = ^dk]. 
In the version published in European Journal of Mechanics, there is an error 
on the equation f3l\) where h{r = 0) = | This is corrected in this postprint. 

When Ni and N2 are calculated using Eqs. ([7]) and (l8|), we find that the 
predicted h{^) variation perfectly agrees with experiments for both directions 
of propagation when ^ < 2.5 (i.e. in the linear regime), as shown in Fig. 4 in the 
particular case uq = 5.3/im, T — Tc = 3.5K, i.e. Bo = 0.015 (dashed straight 
line). As shown in Fig. 4, the numerical predictions h{^) for both upward 
(dot-dashed line) and downward propagations (solid line) are also found to 
agree with both the linear theory and experiments. In addition, Fig. 5 shows 
good agreement between the experimental interface shape (dark symbols) and 
its prediction using both the linear theory (open square) and the numerical 
simulation (solid line) at ^ = 1.35 for both the upward (Fig. 5a) and the 
downward propagation (Fig. 5b). This validates the numerical code versus 
linear analytical predictions as well as its accuracy versus experimental results. 



4.2 Steady-state deformations, nonlinear regime 



Experimentally, for ^ > 2.5, h{^) gradually deviates from linearity. For an 
upward propagation (open circles in Fig. 4) h{^) increases with a larger slope 
than that predicted by the linear theory when 2.5 < ^ < 3.4. For larger values 
of ^, h{^) increases again nonlinearly with a slope lying between that of the 
linear regime and the above mentioned nonlinear one. 
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Conversely, in the case of downward propagation (dark squares), h increases 
with ^ up to ^ ~ 3.4, where its slope diverges. The divergence of h'{^) leads to 
the jetting instability studied in details in [13|; l32l l and illustrated in Fig. 
1 . 

The invariance breakdown of the hump height below the instability threshold 
can be qualitatively explained by considering how the electromagnetic pres- 



sure n(r, 9i, Ot) varies with the incidence angle 6i for both propagations |12 
Whereas the electromagnetic pressure, which varies as ^^fj^, monotonously de- 
creases with 9i in the case of upward propagation because 6t = arcsin(^ sin6'j), 
conversely it continuously increases with 6i in the case of downward propa- 
gation because 9t = arcsin(^ sin 6'j), as shown in the inset in Fig. 6. Thus, 
a downward propagating laser beam exerts an electromagnetic pressure on a 
strongly deformed interface (along which sin 6*4 can reach values comparable 
to unity) larger than an upward propagating beam. The jetting instability 
occurring for the downward propagation can be explained by the fact that the 
beam propagates in this case from the large to the low refractive fluid which 
makes total reflection of light at the interface achievable 



13; 14 



The comparison between numerical simulation results and experimental data 
in the nonlinear regime indicates that, in the case of the upward propagation, 
a satisfactory agreement is found regarding h{^) (see Fig. 4) for the exper- 
imentally investigated values of ^. Moreover, good agreement regarding the 
interface shape is observed up to ~ 5, as shown in Fig. 5a for the particular 
value of ^ = 2.9 and ^ = 4.6. This represents a noticeable progress, since previ- 



ous work on the subject [22|| qualitatively reproduced interface shapes (in the 
case of upward propagation) but did not compare experiments to numerical 
predictions. 

For ^ larger than 3.5, nipple-like interface shapes, shown in Fig. 1, are ex- 
perimentally observed. They are not reproduced by the numerical simulations 
suggesting an additionnal coupling or feedback effect between the exciting 
beam and the soft interface that is not investigated in the present work. As 
also shown in Fig. 4, in the case of downward propagation, the numerical 
simulation qualitatively reproduces the experimentally observed monotonous 
behavior of h{^) up to the jetting instability occurring at ^ = ^f""^, as well as 
the divergence of its slope at the instability threshold. Instability of the inter- 
face is numerically predicted to occur when the wave undergoes total reflection 
and is a priori focused toward the hump tip, i.e. when 9i reaches 9tr at the 
inflection point of the interface. As a matter of fact, since total reflection and 
focusing of the incident wave by the interface is assumed to be responsible 



for the jetting instability [l2||, we assume that beyond the interface will 
actually become unstable. Moreover, beyond ^^""^ the physical model used in 
the numerical model fails, in particular due to Eq. (l23l) because the numerical 
simulation of the jet formation and stability should include the description of 
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the complex interplay between wave propagation and interface deformation 
acting as a wave guide, which is beyond the scope of the present investigation. 
Still, the discrepancy between the values of ^^""^ at which total reflection oc- 
curs along the interface (^^""^ = 3.9 in Fig. 4) and ^1^^ just beyond which the 
bell-shaped interface actually loses its stability (^^ = 4.43 in Fig. 4) calls for 
a study of the interface shape. 

Upward propagation 



N -3 



1 1 1 1 1 1 1 1 1 


' 1 ' 1 ' 

^=1.35 ^jocffl* 
> ^=2 8 

V" '■/ 

\* V 
W '7 

IB J 


' 1 ' 1 ' 1 


.; 

X ■ 




- (a) 




w 




■ Exp 

- Num 

a Linear model 





_7 _6 -5 -4 -3 -2 -1 1 2 3 4 5 6 7 



r 



Downward propagation 




□ Linear model 



-7 -6 -5 -4 -3 -2 -1 1 2 3 4 5 6 7 

r 

Fig. 5. Comparison between experimental and numerical interface profiles for upward 
(a) and downward (b) propagations, = 5.3/im and T — Tc = 3.5if {Bo = 0.015). 
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Figure 5b shows that there is also a noticeable discrepancy between the exper- 
imental and numerically predicted interface shapes for the case of the down- 
ward propagation in the nonlinear regime. More precisely, the experimental 
hump is wider than the numerically predicted one, and the experimental slope 
of the interface shape is smaller than the numerically predicted one all along 
the interface, even when the experimental hump height is well predicted nu- 
merically, as it is the case for ^ = 3.75 and shown in Fig. 5b. This observation 
can explain the discrepancy observed on /i(^) obtained experimentally and 
predicted numerically. As a matter of fact, since the actual slope of the inter- 
face shape is all along the interface smaller than numerically predicted, the 
actual incidence angle 9i is also smaller than that numerically predicted. Con- 
sequently, the total reflection of the wave by the interface reaching a slope 
equal to 9tr and the associated jetting instability are expected to experimen- 
tally occur at a value of ^ larger than that numerically predicted. 



Nevertheless, the reason for this discrepancy between the actual interface 
shape and its numerical prediction in the case of a downward propagation 



is still unexplained. In Ref. [l3| several possible causes were discussed and 



discarded. Among them are the effect of thermocapillary flows due to light 
absorption, the role of optical nonlinearities, and the additional electromag- 
netic pressure applied on the hump tip due to the partial reflection of light 
by the interface below the instability threshold. A more probable cause may 
be the viscous stress applied on the interface by the flow induced by the scat- 
tering of light as a result of refractive index inhomogeneities occurring in the 
bulk of each phase [s^]. Indeed, the resulting scattering force is known (i) to 
be oriented in the direction of propagation and (ii) of radial extension much 
larger than the beam diameter. However, the exact effect of such a flow on 
flnite size humps, in particular the interplay between viscous effects on the 
deformed interface and the feed back of the hump shape on the flow field in 
its vicinity, remains misunderstood and calls for further study. 



4-3 Dynamics of hump formation, nonlinear regime 



Previous work has recently shown that in the linear regime the interface dy- 
namics is accurately described by a linear theory of overdamped interfacial 



circular waves [3j|, and that the interface dynamics is independent of the di- 
rection of propagation. 

However, in the nonlinear regime, the interface dynamics actually depends on 
the wave propagation. 



In Fig. 6, the time evolution of the hump height h(t) is plotted for both the 
upward (dot-dashed line) and downward propagations (solid line) for ^ = 
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Fig. 6. Time evolution of the interface hump height for downward and upward prop- 
agation. ujQ = b.Sfim , T — Tc = 3.5K {Bo = 0.015) and ^ = 3.75. Time t is 
dimensional. The inset shows the variation of the radiation pressure normalized by 
its value at normal incidence, versus the incidence angle for both propagations. The 
dashed line (inset) shows the total reflection threshold Otr- 



3.75, a value belonging to the nonlinear regime. We note that (i) both curves 
coincide during the first 500 ms of the dynamics and eventually separates, 
(ii) the transient is significantly shorter in the case of upward propagation 
and the equilibrium state is reached much faster. The first observation can be 
explained by noting that for t < 500 ms both interface deformations are small 
enough to be described by a linear theory which is precisely independent of 
the beam propagation. The second observation can be qualitatively explained 
by considering the increase of the electromagnetic pressure with the incidence 
angle 9i in the case of downward propagation (as shown in the inset of Fig. 6) . 
As 9i increases along the interface the associated overpressure applying on the 
growing hump increases its final height and thus increases the delay for the 
hump to reach equilibrium. On the contrary, the fact that the electromagnetic 
pressure decreases when 9i increases in the case of an upward propagation 
leads to the opposite conclusion. Consequently, the analysis of the dynamics 
shows that, beyond the regime of linear deformations, our numerical approach 
becomes predictive for the complex nonlinear regime, thus calling for new 
experimental investigation and providing indications on how to perform these 
experiments. 
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5 Conclusions 



A Boundary Integral Element Method has been used to simulate the defor- 
mation of a fluid-fluid interface induced by a focused laser beam. Since inter- 
face deformations with azimuthal invariance were expected, an axisymmetric 
model has been developed. Comparisons between numerical predictions and 
experimental data showed two main results. First, our model can satisfac- 
torily reproduce the variations of interface equilibrium hump height h with 
respect to the ratio of electromagnetic to capillary pressure ^ in regimes of 
small as well as of large deformations. In the case of an upward propagation, 
nipple-like interface shapes were not reproduced, while in the case of a down- 
ward propagation, instability threshold was overestimated probably because 
actual interface shapes are wider than numerically predicted. This could be 
due to bulk steady flows induced by light scattering as a result of liquid index 
inhomogeneitiy which deserves further experimental characterization before 
numerical implementation. A study of the dynamics of the interface for both 
directions of propagation was also performed showing a shorter characteris- 
tic time for the interface reaching its steady state in the case of an upward 
propagation. Future work should take into account the total reflection and 
focusing of light by the deformed interface in order to potentially predict hy- 
drodynamic instability leading to jet formation and micro-droplet emission at 
the tip. Another interesting investigation would be to study the effect of light 



scattering induced flow [33|| on jetting instability. Our flrst results are thus 
very encouraging since they illustrate the accuracy of our numerical model to 
describe and predict intriguing properties of nonlinear behaviors of the cou- 
pling between light and liquid interfaces, a subject of increasing interest due 
to its wide range of applications in soft matter physics on the one hand and 
scarce theoretical results on the other hand. 
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